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^ ■ Abstract 
H 

A new procedure is constructed by means of APS in APLAN language. The procedure solves 
the initial-value problem for linear differential equations of order k with polynomial coefficients 
and regular singularity in the initialization point in the interval [a, 6] and computes the algebraic 
. polynomial y n of given order n. A new algorithm of Lanczos r-method is built for this procedure, 

the solution existence y n of the initial- value problem proved on this algorithm and also is proved 
the optimality by precision of order k derivative of the initial- value problem solution. 
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1. Introduction 



' Problem. To build the procedure with the following parameters: 

. Input: 1. LDUMK — Linear differential equation of order k 

LDUMK := (D[y] = 0); D[y] = A * yW + ■ ■ ■ + C * y + G; (1) 

O ' 

The coefficients A, . . . ,C , G of this equation are known polynomials of independent variable x. 
The solution y of this equation is a function of variable x. The zero point is a regular special point 
^ | of equation (1): 

— the polynomial A becomes zero in a zero point 



A(0) = 0, (2) 



^ ' — the equation solution y is a function analytical in the zero point. 

2. Initial conditions in the zero point 

init_cond(y,d) := {y{d) = Y , y'{d) =Y 1 , y®(d) = Y t }; d := 0; 

These conditions are given in a form of Taylor series coefficients partial sum in the zero point 

T = T,[y] = 2/(0) + y'(0) * x + • • • + y (/) (0)//! * x l - (3) 

The initial- value problem for equations (1), (2), (3) meets the following conditions: 

— the problem has a unique solution 

— the equation 

D[T + V s [u}} =0, (4) 



where 



X 



V[u] := / u(t)dt; 
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s := max{l + 1, k}; (5) 

I ■= deg(T); (6) 

k : = ord-equ(D[y]); (7) 

— equivalent to the problem — has a recurrence relation system for defining the Taylor coefficients 
for the function 

u = c + c\ * x H h c n * x n H 

of the form 

{ coe/Tm/Z(£>[T + y s [u]],i) =0, % = 0, . . . , m, . . . } := Si&S 2 ; 

relating to the function u := y^; 

where Si — a system of r equations like = 0; 

S2 — a system of equations without equations like = 0. 

If the initial-value problem for equations (1), (2), (3) meets these conditions, the solution of this 
problem is 

y := solve(LDUMK, T) = T + V s [u}; (8) 

3. The approximation interval [a, b}. Initial point d € [a, b]. 

4. The order n of initial-value problem (1), (2), (3) solution (8) approximation sought within 
the interval [a, b] is the polynomial 

y n = opt. T — met(LDUMK, T,[a,b],n) = c + c\ * x H \- c n * x n (9) 

Output. The polynomial (9). This polynomial meets the initial conditions (3) 

y n = T + V s [u p }; (10) 

where 

u P = yi s) ; p:=n-s; (11) 

— a polynomial that approximates the equation (4) solution u = y^ . 

The polynomial y n (9) derivative of order k within the interval [a, b] approximates the initial- 
value problem (1), (2), (3) solution (8) derivative of order k optimally by precision 



\y {k) -yW\\ L2{a , b ; p) /E[n-k,y( k \L 2 (a, b; p)} < Const < 00, (12) 



where 



ii2 

\ V \\L 2 (a,b;p) 



f v\x) p(z(x)) dx; p(z) := (1 - z 2 )^ - 1/2); (13) 

J a 

2!(x):=2*(x-o)/(6-o)-l; (14) 
E[n-k,y^,H]:= inf \\yW _ + . . . + Cji _ k * x n - k )\\ H ; (15) 

- the value y^ of best function approximation by the algebraic polynomials of order n — k in the 
space H 

If the polynomial y n meets the inequation (12) and the initial conditions (3), (10), it is optimal 
by order the joint problem (1), (2), (3) solution y (8) approximation and its derivatives of order 
i = 1, 2, . . . , k — the following identities take place 



y 



®-V$ = V k - i [v-v n - k ]; i = 0,1,2,..., k 
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where v := y {k >; v n - k ■= Vn and ||V||q 0]i)] = max{|a|, \b\}. 

Relevance of the problem. The equations (1), (2) define the most part of special functions of 
mathematical physics. The order I (6) of polynomial T (3) — the conditions of initial- value problem 
for equations (1), (2) — usually is not equal to k — 1, where k — the order (7) of equation (1). The 
Bessel function of order is defined uniquely by the linear differential equation with polynomial 
coefficients of order 2 and the polynomial (3) of order 0. The Frenel integrals are defined uniquely 
by the linear differential equation with polynomial coefficients of order 3 and the polynomial (3) of 
order 3. 

Linear differential equations with polynomial coefficients comprise the mathematical modelling 
apparatus for physical and technical processes [1]. The computer algebra systems already became 
the natural mathematical modelling media. While analyzing the processes described by equation (1) 
the computer algebra systems usually perform analytical transformations of equation (1) solution. 
The equation (1) solution usually is not a composition of functions, that can be symbolically 
transformed by the computer algebra system. Thus the computer algebra system transforms the 
equation (1) solution approximation. Usually this approximation is a polynomial. The computer 
algebra systems have the efficient programming means of symbolic polynomial transformations. 

The main criterion for modeling is the model precision. According to this criterion the polyno- 
mial should approximate the solution optimally by precision order. The computer algebra systems 
have the procedure for computing the Taylor series of the initial- value problem solution for ordinary 
differential equations. The Taylor series is not the optimal function approximation apparatus - 
the order of function y Taylor series remainder term is bigger than the value of this function best 
approximation 

\\y-T n [y}\\ C[aM /E[n,y,C [aM }=0(q n ), q>2 

Lanczos [1] developed the r-method of solving the initial- value problem for equation (1). Luke 
applied the Lanczos r-method for computing the Fourier-Chebyshev coefficients for special functions 
of mathematical physics. These coefficients are the foundation for the procedures of computing the 
special functions of mathematical physics in mathematical computer software. The importance of 
these procedures initiated the development of new approximative methods of solving the initial- 
value problem for the equation (1) — Clenslaw method, Miller method, V. K. Dzyadyk a-method 
[2] and others. 

V. K. Dzyadyk a-method [2] solves the the initial- value problem for the equations (1), (3) without 
the singularity (2) optimally by precision order. The analytical methods of solving the initial-value 
problem for equations (1), (2), (3) optimally by precision are not developed yet. 

2. APLAN-procedure for Lanczos r-method 

Data structure. 

Input: 1. The equation (1) is 

LDUMK : = ( A * dif ( y , k ) + . . . + C*y + G = 0); 

where y — an atom, coefficients A , . . . , C, G — polynomials of equation (atom) x. These 
polynomials are separate terms (they are not a multiplication of terms). Usually they have the 
natural form for mathematics and are taken in parenthesis. 
2. Polynomial T (3) has natural form for mathematics 

T:=d+e*x+. . . + f * x ~ q ; 
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The argument x of this polynomial is an atom. 

3. The approximation interval [a, b] description is a list of interval ends. 

interval := (a , b) ; 

Output. The procedure computes the polynomial y n (9) with numeric coefficients. The form 
of this polynomial is identical to the form of polynomial T. 
Algorithm 1. 

1. To compute the equation (1) operator D[y]. 

2. To compute the differential operator D[y] order k (7). 

3. To compute the order I (6) of the polynomial T (3). 

4. To compute the algorithm parameter s (5). 

5. To compute the order p (11) of the polynomial u p . 

6. To compute the polynomial of order p with symbol coefficients 

Up = co + ci * x + • • • + Cp * x p ; (16) 

7. To compute the polynomial u p (16) transformation (10). 

8. To compute the transformation of polynomial y n (16), (10) by operator D[y] (1) 

D[y n ] = D[T + V s [up\], (17) 

9. To compute the zero order of the polynomial D[y n ] (17) in the zero point 

r = deg.nul{D[T + V s [u p }}); (18) 

10. To compute the polynomial D[y n ] (17) regularization 

D [y n ] =D[T+V s [u p ]}/x r ; (19) 

the polynomial Do[y n ] does not have zeroes in zero point. 

11. To compute the polynomial D o[y n ] (19) order 

m = deg(D[T + V s [u p }]/x r ); (20) 

12. To compute the auxiliary V. K. Dzyadyk a-method [2] polynomial of power m (20) and 
order p (11) (the discrepancy) within the interval [—1, 1] 

E m (x) = E(m,p, x) =n* f p +i(x) H h r m _ p * f m {x); (21) 

The discrepancy basis consists of Chebyshev polynomials of first type [2] 

fi(x) = cheb(i,x) = cos{i * arccos(x)); (22) 

13. To compute the linear transformation (14) of the interval [a, b] into the interval [—1,1]. 

14. To compute the auxiliary V. K. Dzyadyk a-method [2] polynomial (21) transfer z(x) (14) 
onto the interval [a, b] 

E m (z) = E(m,p, z) = subs(x = z, E(m,p, x)); z = z(x); (23) 

15. To compute the V. K. Dzyadyk a-method [2] approximation of the left side of the regularized 
equation (4) — the sum of polynomials (19), (23) 

D[T + V s [u p ]]/x r + E m {z) (24) 
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16. To compute the system of linear algebraic equations (SLAE) 



S = { coefTayl(D[T + V s [u p ]]/x r + E m {z),i) = 0, % = 0, . 



,m }; 



(25) 



The system variables are the coefficients of the polynomials u p , E, 



•m 



c , . . . ,c p , 



T~l 5 • • • , T m —p 



(26) 



17. To compute the SLAE (25) solution — the coefficients (26) values 



Coef := solve(S) := { c = d, . . . , c p = e, n = /,... , r m _ p 



9h 



(27) 



18. To compute the polynomial u p (11) with numeric coefficients 



u, 



v 



ser{ Coef, p ) := d + ■ ■ ■ + e * x p ; 



(28) 



The values of the polynomial u p coefficients (26) are defined by the identities (27). 

19. To compute the polynomial (28) transformation (10). This polynomial is the initial-value 
problem (1), (2), (3) solution (8) approximation sought within the interval [a, b}. 

Algebraic specification of algorithm 1. 



u_p := ser(p , Coef); /* approximation of dif(y,s) */ 
y_n := T + n_int( u_p , s ); /* approximation of y */ 

The structure of computations results obtained by procedure operators. 

The procedure solves SLAE (25) relative to the coefficients (26) of form 

c(0) , . . . , c(p) , c(p + 1) , . . . , c(m) 

These coefficients are indexed atoms. The procedure transforms the polynomials u p (16) E m (21) 
with these symbol coefficients to obtain the system (25). 
The polynomial (16) with symbol coefficients looks like 

u_p :=c(0)+c(l)*x+. . .+c(p)*x"p; 



let( LDUMK , Dy = ) ; /* operator Dy */ 

k := ord_equ( Dy ) ; /* order Dy */ 

s := deg( T ) + 1; /* order T + 1 */ 

(k>s)->(s:=k); /* parameter of the method - s */ 

p := n - s; /* power of the polynomial u_p */ 

u_p := main_pol(p); /* u_p with coefficients c(i) */ 

y_n := T + n_int( u_p , s ); /* T + V~s[u_p] */ 

Dn := canplf (sub_du(Dy , y_n)); /* polynomial D[y_n] */ 

r := deg_nul (canplf (ein_pol (Dn) )) ; /* zeroes of D[y_n] */ 

Dn — > ndiv_x ( Dn , r ) ; /* Do[y_n] */ 

m := deg(canplf (ein_pol(Dn) ) ) ; /* order Do [y_n] */ 

Em := Enl(p, m-p) ; /* E_m for [-1,1] */ 

z — > canplf ( -1 + (2/(arg(interval,2) + (-1) * 

arg (interval , 1) ) * (x + (-1) * arg (interval , 1) ) ); 

Em — > canplf ( subs( x = z, Em )); /* E_m(z) */ 

Dn — > canplf (Dn + Em) ; /* Do [y_n] + E_m */ 

S := pol_equ(Dn , m) ; /* SLAE - problem approximation */ 

Xn := c; Coef := solve (S); /* SLAE solution */ 
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The result of transforming this polynomial by operators 

T + n-int(ujp, s), canplf (sub-du(Dy,yjn)), ndiv-x(Dn,r) 

is a sum of addends like c(i) * x^$b , where $ — the APLAN operation of multiplying the term by 
constant, constant b is rational and looks like rat(p,q) -=p/q ', 
The discrepancy (21) with symbol coefficients looks like 

E_m := c(p + 1) * cheb(p + l,x) + . . . + c(m) * cheb(m,x) ; 

The result of transforming this polynomial by operator 
canplf (subs(x = z,Em)) is a sum of addends like c(i) * x 3 $b . 
Thus SLAE (25) is a list of equations of form 

S :=(... , c(m) $ f + . . . + c(0) $e+d=0 , ... ); 

The solution of SLAE (25) by the procedure is a list of identities 

Coef := ( c(0) = d , ... , c(p) = f , ... , c(m) = g ); 

The procedure transforms the list of identities Coef into the polynomial u p with numeric natural 
form coefficients 

u_p := d + e * x + . . . + f * x ~ p ; 

Conclusions from the algebraic specification of the procedure and data structure: 

1. The procedure has known operators [3]. These operators perform computations in rational 
numbers arithmetics. The length of numerator and denominator of these numbers is not limited. 
Thus the procedure operators do not add into the computations result the errors of performing the 
floating point operations. 

2. The procedure has the polynomial complexity by parameter n 

0{n*Q(canplf,n 2 )) + 0(n 3 ), (29) 

where Q(canplf, m) — the complexity of transforming by operator canplf the polynomial P, m — 
number of addends of the polynomial canplf (P). The operator canplf is the internal operator of 
the APS solver gr_solve.exe. It reduces the polynomials into the canonical form for the APS. If P 
is polynomial transformed by the procedure, then canplf (P) is a sum of addends like c(i) * x J $6 . 
Computational experiment with the procedure. 
The initial-value problem 

x*y" -y' + 4*x 3 *y = 0, T = x 2 ; (30) 

has the unique solution y = sin(x) 2 . The APLAN-description of initial-value problem (30) and the 
interval [—1,1] for the procedure implemented above is 

process [1] := ( LDUMK : = 

( x * dif(y , 2) + (-1) * dif(y , 1) + (4 * x~3) * y = ); 
T := x " 2 ; 
interval := (-1 , 1) ; ...); 

The results of this initial- value problem transformation by the procedure built with the param- 
eter value n = 4 are: 
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let( LDUMK , Dy = ) ; Dy := x * dif (y , 2) + 

(-1) * dif(y , 1) + (4 * x " 3) * y ; 
k := ord_equ( Dy ) = 2 ; 
s := deg(T) +1=3; 
(k>s)->(s:=k); s:=3; 
p : = n - s = 1 ; 

u_p := main_pol( p) =cO+cl*x; 
y_n := T + n_int( u_p , s ) = 

c 1 * x " 4 $ rat (1,24) + c0*x~3$ rat (1,6) + x " 2 ; 
Dn := canplf ( sub_du( Dy , y_n ) ) = 

c 1 * x ~ 7 $ rat (1,6) + cl*x~3$ rat (1,3) + 
c * x~6 $ rat (2 ,3) + c * x~2 $ rat (1,2) + x~5 $ 4 ; 
r := deg_nul( Dn ) := 2 ; 
Dn — > ndiv_x ( Dn , r ) = 

c 1 * x ~ 5 $ rat (1,6) + c 1 * x $ rat (1,3) + 
c * x ~ 4 $ rat (2, 3) +x~3$4+c0$ rat (1,2); 
m := deg( canplf ( ein_pol( Dn ) ) ) = 5 ; 
Em := Enl( p , m - p ) = 

c 5 * x A 5 $ 16 + c 5 * x * 3 $ -20 + c5*x$5 + 

c4*x~4$8+c4*x~2$-8+c3*x~3$4 

+c3*x$-3+c2 *x~2$2+c4+c2$-l; 
z := x ; 

Em — > canplf ( subs( x = z , Em ) ) := 

c 5 * x A 5 $ 16 + c 5 * x * 3 $ -20 +c5*x$5+ 

c4*x~4$8+c4*x~2$-8+c3*x~3$4 

+c3*x$-3+c2 *x~2$2+c4+c2$-l; 
Dn — > canplf ( Dn + Em ) = 

c 5 * x A 5 $ 16 + c 5 * x * 3 $ -20 +c5*x$5+ 

c4*x~4$8+c4*x~2$-8+c3*x~3$4 

+c3*x$-3+c2 *x"2$2+c4+c2$-l+ 

c 1 * x " 5 $ rat (1,6) + c 1 * x $ rat(l,3) + 

c * x " 4 $ rat (2, 3) +x~3$4+c0$ rat (1,2); 
S := pol_equ(Dn , m) = ( 

c4+c2$-l+c0$ rat (1,2) = , 

c5$5+c3$-3+cl$ rat(l,3) = , 

c4$-8+c2$2=0 , 

c 5 $ -20 +c3$4+4=0 , 

c4$8+c0$ rat (2, 3) = , 

c 5 $ 16 + c 1 $ rat (1,6) = ) ; 
Coef := solve (S) =(c4=0,c5= rat (1,14) , c 2 = , 

c 3 = rat (-9, 14) , c = , c 1 = rat (-48, 7) ) ; 
u_p := ser(p , Coef) = rat(-48 , 7) * x ; 

y_n := T + n_int( u_p ,s)=x~2+x~4$ rat(-l , 7); 

3. The solution existence by the procedure 

The Hilbert space H. The equation (1) solution is the function analytical throughout the 
whole complex plane. It has the finite number of poles. Only the zeroes of the LDEPC (1) 
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coefficients could be the poles of this function. Thus it is natural to consider the equation (4) and 
its approximation (25) in the space of functions analytical in the zero point. We will consider the 
algorithm 1, where the basis meets the following conditions: 

— the basis elements are the polynomials of power i 

fo(x), /i(x), / 2 (x),... deg(fi(x)) = i; (31) 

- the analytical functions have the Fourier series on the basis (31) within the interval [a, b] 

u = c * fo(z(x)) H hc m * f m (z(x)) H ; 

v = d * fo(z(x)) H h d m * f m {z{x)) H ; 

Thus basis (31) defines the Hilbert space H in the set of analytical functions. This space has the 
scalar product and norm 

(u, v ) H = c * d H \-c m *d m -\ \\u\\ 2 H = (u,u) H (32) 

The basis (22) of V. K. Dzyadyk a-method is is a basis (31) special case. It defines the Hilbert 
space L 2 (— l,l;p). Linear transfer (14) of this basis onto the the interval [a,b] — {f m (z(x))} - 
defines the space L 2 (o, b; p). 

The space of algebraic polynomials of order p 

H p := S P [H}; S p [c * f H h c p * f p H ] := c * fo H h c p * f p ; 

— projecting operator by basis (31), is a subspace of the space H. 

Theorem 1. Let SLAE (25) have the discrepancy basis (31), p > po and equation (4) linear 
regularization operator 

(D[T + V s [u]]/x r = 0) :=(£[«] + / = 0); (33) 

in the Hilbert space H (32) meet the following conditions: 

— operator L definitional domain D(L) is dense in the space H, 

— operator L domain of values R(L) is dense in the space H, 

- the operator L transforms the domain D(L) into R(L) biuniquely, 

- the subspaces H p — of algebraic polynomials of order p — H p and L[H p ] are closed in the 
H space, 

— the subspaces sequence L[H p ] is dense in H down to the limit; 

lim 7(p) >0, 7(p) := inf \\S p [z p ]\\ H ; 

P-^oo z p eL[H p ], \\z p \\ H =l 

then the SLAE (25) solution Coef (27) exists and it is unique. 

Proof layout. We proved the equivalence of the SLAE (25) with the basis (31) with equation 
(33) approximation by projective method with projecting operator by basis (31). Thus the theorem 
1 is the consequence of main theorem of projectional method convergence. 

4. The solution optimality by the procedure 

To prove the inequation (12) we transform the integral equation (4) into the linear integral equation 
of the third type 

subs(y = T k _ 1 [y] + V k [v],D[y] = 0) := {A * v + • • • + C * V k [v] + g = 0); (34) 
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The equation (34) solution is the function 

v = y^ :=V s - k [u}+T^; 
The equation (34) is equivalent to the linear integral equation of the second type 

M[v}+ g/A = 0, M[v] := v + ■ ■ ■ + C /A * V k [v}; (35) 
and the initial conditions 

T s ^M=T {k) 

Similar transformation of the integral equation (4) approximation (24) by V. K. Dzyadyk a- 
method lead to the equation 

M[v n - k ] + g/A + x r /A*E(m,p,z(x)) = 0; (36) 

where E(m,p,x) (21) is the auxiliary V. K. Dzyadyk a-method [2] polynomial with the basis (31), 
and the initial conditions 

r s _ fe _ 1 K_ fc ] = r( fc ) 

Equation (36) is the equation (35) approximation by addition of the discrepancy x r /A*E(m,p, z[x)). 
This discrepancy basis is the transformation of basis (31). 

/„(*(*)) * x r /A, h{z{x)) * x r /A, f 2 (z(x)) * x r /A, . . . (37) 

Theorem 2. Let: 

- linear operator M[v] of the equation (35), M : H — > H, where H — the Hilbert space 
of analytical functions with the scalar product (32) by the basis (31) with the initial conditions 
T s -k-i[v] = have reverse operator M^ 1 ; 

- the operator M~ l (36) be uniformly limited on the basis (37) functions 

\\M- 1 [f i (z(x))*x r /A]\\ H <W, i = 1,1 + 1,...; (38) 

- SLAE (25) with the basis (31) have the unique solution (27); 
-for n = I + s + 1, 1 + s + 2, . . . the matrices 

{Qn} = {(M-V'M * f p+ i{z(x))],f p+3 (z{x))) H , i,j = l,...,m-p}, (39) 

where p = n — s, m = deg((D[V s [u p ] + T]) — r, r = degjnul(D\V s [u p ] + T\); 

have their reverse {Qn} -1 and the reverse matrices norms be uniformly limited by n 

||{Qn} _1 || := , , max ||{Qn} _1 {Ci, • • • , Cm- P } T \\i 2 < Q (40) 

|Cl|+-+|Cm-p| = l 

Then assuming n > I + s the inequation (12) is valid and its constant is 

Const = Q*W 



Proof layout. The dependence of equation (35) solution v approximation by the equation 
(36) solution v n ^.k from the equation (36) discrepancy coefficients 

\\v -v n -k\\ H < W * (|ri | H h \T m - p \) 
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and the dependence of discrepancy coefficients from Fourier coefficients of the equation (35) solution 
v by the discrepancy basis 

(|n(p)| H h |r m _ p (p)|) 2 < Q 2 * ((v,f n _ k+1 ) 2 H H h {v,f n ^k+m-p) 2 H) 

are established. 

Thus the inequation (12) appears directly from these inequations. 

Example 1. We solved the initial- value problem by the procedure built within the interval 
[—1,1] with n = 4 , 6 , . . . , 22 and computed the polynomials 

Vn = opt. T — method(x * y" — y + 4 * x 3 * y = 0, y(0) = 0, T = x 2 , [—1, 1], n); 

For these polynomials the equation (12) constant Const estimation in the space C = (^[-1,1] has 
the following values 

{ \\sin"(x 2 )-y'; i \\c/E[n-2, S in"(x 2 ),C], n = 4, 6, . . . , 22; } = 

{ 1.6 ,1.5 ,2.1 ,1.5 ,2 ,1.7 ,2 ,1.7 ,2 ,1.8 ,1.9 }; 

5. Postamble 

A new algorithm of solving the initial- value problem for the equations (1), (2), (3) by Lanczos r- 
method is built in the article. The APLAN language procedure is built according to this algorithm 
by the programming technology in the APS. This procedure solves the initial- value problem for 
the equations (1), (2), (3) optimally by precision order. It illustrates the high efficiency of APS 
toolkit [3] for creating the new optimal analytical approximative methods of solving the functional 
equations. 
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